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1  Introduction 

The  goals  of  this  research  is  to  overcome  the  difficulties  created  by  mesh  shearing 
in  the  simulation  of  aircraft  with  large  control  surface  deflection.  The  origin  of 
the  difficulties  for  the  correct  modeling  of  large  control  surface  deflection  is 
that  they  introduce  shearing  in  the  mesh  at  the  place  where  the  control  surface 
meets  the  fixed  part  of  the  wing  (see  Fig.  1).  Such  a  geometrical  shearing  is  very 
difficult  to  handle,  in  particular  in  Finite  Volume  or  Finite  Element  approaches. 
One  option  to  tackle  this  problem  is  the  use  of  Chimera  grids. However  Chimera 
grid  approaches  are  relatively  expensive  in  geometric  computation  and  introduce 
an  interpolation  error  due  to  the  overlap  of  the  grids.  We  have  proposed  to 
explore  two  alternate  methods  to  tackle  this  issue: 

1.  Meshfree  Methods. 

2.  Level-Set  Methods. 

During  the  first  year  of  this  research,  we  developed  the  Meshfree  method  for 
aeroelastic  problems  and  created  a  new  approach  to  dealing  with  Neumann 
Boundary  Conditions.  We  then  proceeded  to  implement  this  approach  in  two 
dimensions  and  showed  initial  aeroelastic  results  for  an  airfoil  moving  in  a  forced 
pitching  motion.  The  experiments  showed  that  the  method  worked,  though  the 
cost  of  the  method  tended  to  indicate  that  its  use  should  be  limited  to  areas 
where  other  methods  would  not  be  appropriate. 

In  the  second  year,  we  then  tackled  the  problem  of  Control  Surface  Deploy¬ 
ment.  This  efTort  resulted  in  some  partial  success.  However,  new  problems 
surfaced  with  conditioning  of  the  resulting  system  and  we  started  the  level-set. 
type  approach  to  avoid  the  ill-conditioning  problem.  This  approach  has  shown 
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Figure  1:  Mesh  Shearing  Regions  due  to  Control  Surface  Deflections 


good  results  for  any  geometry  and  at  the  end  of  the  project,  work  was  continu¬ 
ing  in  implementing  the  full  time  dependent  solution  of  problems  with  moving 
control  surfaces. 


2  The  Meshless  Approach  to  Aeroelasticity 

2.1  The  reasons  for  a  Meshless  Approach  to  Aeroelasticity 

We  have  chosen  to  approach  the  problem  of  simulating  aircraft  with  large  con¬ 
trol  surface  deflection  with  a  Mcshlcss  Approach.  Mcshless  methods  have  been 
applied  to  problems  from  Crack  propagation  to  Astrophysics  and  Fluid  prob¬ 
lems.  In  general,  such  methods  are  particularly  adapted  where  the  boundary  of 
the  computational  domain  greatly  evolves  during  the  time  of  the  simulation. 

Meshless  Methods  are  particularly  adapted  to  problems  where  the  boundary 
significantly  evolves  in  time,  because  they  never  require  to  form  or  maintain  a 
mesh  topology.  Unlike  Finite  Element  or  Finite  Volume  methods,  where  the 
motion  of  nodes  is  constrained  by  the  necessity  to  preserve  positivity  of  the 
volumes  of  the  discretization  elements,  the  Meshless  Methods  can  accept  almost 
arbitrary  relative  motions  of  all  the  nodes. 

The  only  constrain  imposed  on  the  relative  motion  of  the  nodes  is  that  a 
sufficient  number  of  nodes  are  close  to  any  point  in  the  domains,  so  that  the 
interpolation  functions  are  well  defined  for  every  point  in  the  domain.  The  exact 
number  of  nodes  and  the  definition  of  close  depends  on  the  particular  variant 
of  Meshless  Methods  being  used.  In  general,  such  constraints  are  rather  weak 
and  can  easily  be  enforced  with  no  significant  difficulty. 

These  facts  are  major  advantages  for  the  Meshless  Approaches  when  applied 
to  Aeroelasticity  problems.  This  is  true,  not  only  when  large  control  surface 
deflections  arc  present.,  but  also  when  large  relative  motions  of  objects  have 
to  be  simulated.  Examples  of  such  problems  are  Missile  Separation  problems, 
Ejection  Seat  simulations  and  Formation  Flying  of  several  aircraft. 
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2.2  Advantages  of  Meshless  Approaches  to  Aeroelasticity 

•  Ease  of  handling  arbitrary  relative  motion  of  boundaries.  Because  mesh- 
less  methods  do  not  require  to  maintain  a  topological  coherence  to  the 
mesh,  nodes  can  be  moved  almost  arbitrarily  relative  to  one  another. 

•  Ease  of  refining  the  solution  by  inserting  new  nodes.  New  nodes  can 
easily  inserted  to  refine  the  solution  where  high  gradients  are  observed. 
The  insertion  does  not  require,  as  in  FE  or  FV  methods  to  construct  a 
new  mesh  topology. 

2.3  Disadvantages  of  Meshless  Approaches  to  Aeroelastic¬ 
ity 

The  major  disadvantages  of  Meshless  Approaches  when  applied  to  Aeroelastic 
problems  are: 

•  Cost:  The  computation  of  the  shape  functions  with  some  methods  (for  ex¬ 
ample  RKPM)  is  much  more  expensive  than  a  Finite  Element  equivalent. 

•  Essential  Boundary  Conditions:  Imposing  Essential  Boundary  Conditions 
in  Meshless  Methods  is  not  a  trivial  matter.  This  is  shared  with  all  Mesh¬ 
less  methods  and  is  due  to  the  fact  that  the  value  of  the  solution  at  a 
computational  node  is  not  equal  to  the  value  of  the  node.  Or  phrased 
mathematically,  the  Meshless  shape  functions  do  not  exhibit  a  Kronecker 
delta  property. 

•  Complexity  of  integration,  choosing  Gauss  Points  and  finding  nodes  with 
influence  over  the  Gauss  Points. 

3  Results  on  Meshless  Methods 

A  graduate  student,  Vivek  Kaila,  was  recruited  to  work  on  this  project  and 
began  working  March  1st  2004.  Our  approach  is  based  on  the  use  of  the  Repro¬ 
ducing  Kernel  Particle  Method  combined  with  the  Streamline  Upwind  Petrov- 
Galerkin  (SUPG)  method. 

3.1  RKPM  Approach 

Of  the  various  methods  that  exist  in  the  literature  on  Meshless  Methods,  we 
chose  to  concentrate  on  the  Reproducing  Kernel  Particle  Method  (RKPM),  due 
its  use  by  other  researchers  in  CFD  problems.  The  research  did  not  seek  to  be 
on  the  chosen  Meshless  Method,  but  on  the  adaptation  of  Meshless  Methods  to 
Aeroelastic  problems.  The  developments  made  in  this  research  are,  and  should 
remain,  applicable  to  any  Meshless  Method. 
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Common  to  all  Meshless  Methods,  the  solution  is  constructed  from  some 
nodal  values 1  Ui  and  shape  functions  $j(x)  are  constructed  for  each  node.  The 
solution  is  then  given  by: 


(1) 


/ 


The  difference  between  the  different  existing  methods  is  in  the  construction 
and  resulting  properties  of  the  shape  function.  We  note  that  the  difficulty 
in  the  imposition  of  essential  boundary  conditions  is  a  result  of  the  fact  that 
^  6jj  where  xj  is  the  location  of  node  J. 

3.2  Imposing  Boundary  Conditions 

As  was  mentioned  before,  the  main  difficulty  with  Meshless  Methods  is  to  impose 
Essential  Boundary  Conditions.  Several  existing  approaches  were  considered  but 
were  deemed  unsatisfactory: 

•  Lagrange  Multiplier  Methods:  Such  methods  impose  the  essential  bound¬ 
ary  conditions  in  a  weak  manner.  Either  by  collocation  —  in  which  case 
the  boundary  conditions  are  satisfied  only  at  a  discrete  set  of  points  -  or 
in  an  average  form.  The  resulting  solution  is  equivalent  to  having  solved 
the  problem  around  a  slightly  different  structure  than  the  intended  one. 

•  Blending  of  Meshless  and  Finite  Element  methods:  Such  methods  which 
transition  the  Meshless  solution  to  a  Finite  Element  (FE)  solution  on 
the  boundary  do  not  suffer  of  the  aforementioned  problem.  However,  if 
refinement  of  the  solution  is  needed,  they  require  a  rc-meshing  of  the 
FE  mesh  near  the  boundary  to  capture  the  same  fine  solution  as  in  the 
Meshless  Domain.  Such  re-meshing  is  much  more  expensive  to  perform 
than  the  insertion  of  nodes  in  the  Meshless  approach  and  requires  an 
added  interpolation  of  the  solution  between  meshes,  if  this  operation  is 
performed  during  a  time  dependent  simulation. 

We  decided  to  design  a  new  method  for  the  imposition  of  the  Essential 
Boundary  Conditions,  with  the  following  features: 

•  The  method  should  preserve  the  spatial  resolution  of  the  Meshless  solution 
without  severe  difficulties. 

•  It  must  impose  the  boundary  condition  strongly  on  the  entire  boundary. 

•  It  must  be  easy  to  adapt  for  boundary  surfaces  undergoing  arbitrary  and 
potentially  large  relative  motions. 

Our  efforts  have  resulted  in  the  design  of  a  Projection  Based  Boundary 
Condition  treatment  which  we  will  now  describe  below: 

1  These  values  are  associated  with  the  node  but  are  not  the  value  of  the  solution  at  the 


node. 
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Though  we  have  started  with  the  inviscid  flow  condition,  and  will  only 
present  this  case  in  what  follows,  the  method  extends  fairly  easily  to  viscous 
flows. 

In  inviscid  flows,  the  common  type  of  boundary  condition  is  the  slip  con¬ 
dition,  where  the  flow  must  be  tangent  to  the  airfoil  or  wing  surface.  Such  a 
condition,  in  the  presence  of  a  structural  velocity  at  the  interface  us  stipulates 
that  the  flow  velocity  uj  must  be  such  that: 


up  n  =  us  n  on  I>/s  (2) 

Let  us  remind  the  reader  that  the  solution  quantities  are  given  by  the  conser¬ 
vative  variable  vector  which  in  two  dimensions  is  given  by: 


U  = 


(  P\ 

pu 

pv  I 

pE  / 


(3) 


where  u  v  are  the  components  of  the  velocity  vector  and  E  is  the  total  energy 
density. 

We  then  notice  that  given  a  solution  Um  that  would  be  obtained  by  the 
Meshless  Approach,  we  can  build  a  solution  U  that  satisfies  the  Boundary  Con- 
dition  by  the  use  of  a  project: 


U  =  PUm  +  Qus 


(4) 


where 


and 


(10  0  0\ 

0  1  -  nxnx  —nxny  0  | 

0  - nxny  1  -  nyny  0  I 

0  0  0  1/ 


0  0  \ 

^X^X  I 

^X^y  TlyTiy  I 

0  0  / 


It  is  therefore  possible  to  construct  a  continuous  solution  by  gradually  blending 
the  modified  solution  given  by  equation  4  and  the  meshless  solution  UM. 

In  the  presence  of  a  single  surface,  we  define  a  ramping  function  r(x)  which 
is  such  that  r(x)  =  1  on  IV/s  and  becomes  0  away  from  the  surface.  Then  the 
solution  is  given  by: 


U  =  ((  1  -  r(x))I  +  r(x)P)UM  +  r(x)us  (5) 

The  ramping  function  must  transition  between  0  and  1  around  the  interface  in 
a  distance  compatible  with  the  expected  resolution  of  the  Meshless  Method. 

The  method  can  then  be  extended  to  multiple  surfaces  by  creating  one  ramp¬ 
ing  function  that  satisfy  that  the  ramping  function  r,(x)  for  surface  i  is  0  away 
from  that  surface  and  also  on  any  other  boundary  surface,  while  it  is  1  on  its 
respective  surface. 
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Figure  2:  Pressure  contours  at  M=0.6  on  a  2000  node  Mesh 


4  Difficulties  with  the  Meshless  Method 

Though,  overall,  the  Meshless  Method  seems  to  perform  reasonably  well.  How¬ 
ever,  when  control  surfaces  are  deployed,  as  the  nodes  attached  to  the  various 
components  follow  them  in  their  motion,  nodes  attached  to  one  component  can 
get  extremely  close  to  nodes  associated  with  another  component.  As  a  result, 
two  problems  became  evident  as  more  simulations  were  performed  and  time 
dependent  with  control  surface  deployment  problems  were  explored: 

•  on  coarse  mesh,  the  solution  exhibits  unacceptable  oscillations. 

•  as  the  flap  angle  changes,  node  getting  close  create  a  very  ill-conditioned 
system. 

The  second  problem  is  a  major  one  and  can  easily  be  understood  because  the 
equation  for  the  shape  functions  starts  exhibiting  extremely  strong  local  gradi¬ 
ents  when  nodes  get  too  close  to  one  another.  The  first  problem  is  linked  to 
the  second  and  is  exhibited  in  figure  \  Consequently,  it  was  realized  that  some 
method  would  have  to  be  developed  so  that  solution  from  different  meshes  are 
not  active  in  the  same  area  to  avoid  this  problem.  It  seemed  that  a  level-set  type 
method  used  to  blend  the  solution  from  two  meshes  would  provide  a  reasonably 
inexpensive  solution  to  this  problem. 

5  Levelset  Method 

The  levelset  method  is  used  to  blend  the  solution  of  two  meshes  as  shown  in 
figure  5  For  each  mesh,  we  created  two  levelset  methods,  and  ty,: 

•  3>  is  akin  to  the  Distance  from  the  airfoil 

•  <f>,  =  0  on  Ft ;  =  1  on  far-fiold  i 

•  is  akin  to  the  Distance  from  Non-Physical  Mesh  boundary 
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Figure  3:  Cp  contour  on  an  airfoil  with  a  deflected  control  surface 
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Figure  4:  The  airfoil  and  flap  meshes 
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Figure  5:  Level  set  functions 


Both  functions  are  created  by  solving  a  simple  Laplace  problem.  Figure  fig:phipsi 
shows  the  contours  of  such  functions.  The  level-set  functions  are  then  used  to 
blend  the  solutions  U\  and  U2  from  both  meshes:  Let: 


D 1 


-  d  t  1 
c  +  2 


(6) 


and 

D2  =  1  -  A  (7) 

In  these  equations,  d  is  a  function  of  the  angular  deflection  of  the  flap  and 
measures  the  difference  'Ll  —  ty2  at  the  intersection  points  of  the  meshes,  c 
is  a  parameter  to  adjust  the  width  of  the  transition  and  the  function  D}  and 
consequently  D2  are  then  clamped  between  0  and  1  / 

We  then  define  thee  influence  coefficients  for  each  mesh,  s\  and  s2 : 


c  _  _ *2(1  -  $i)/£>2 _  o  ,  _  c 

1  $2(i -*i)/e2  +  *i(i -*2)/z?i  ’  1 

Finally,  we  use  these  functions  to  blend  the  the  two  mesh  solutions: 


U  —  S\U\  +  S2U2 

Figure  5  and  5  shows  the  contours  of  Si  =  0.9,  S\  =  0.5  and  S\  =0.1  from  left 
to  right. 


5.1  Node  activation  and  de-activation 

One  notes  that  for  some  nodes,  the  52  function  associated  with  their  mesh,  will 
be  identically  zero  over  the  whole  domain  of  influence  of  the  node.  Such  nodes 
will  thus  have  no  contribution  to  the  solution  and  are  to  be  de-activated  from 
the  solution  process. 
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Figure  6:  Contours  of  the  S\  function 
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Figure  7:  Close-up  of  the  contours  of  the  S\  function 
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Active  zone 


Inactive  zone 


Space-Time  slab 


Figure  8:  One  dimensional  Space-Time  slabs  associated  with  a  single  mesh 


We  also  note  that,  as  the  angle  of  deflection  of  the  control  surface  varies  in 
time,  nodes  can  become  inactive  or  active  over  time.  It  is  therefore  necessary 
to  devise  a  procedure  to  deal  with  this  issue.  The  approach  we  devised  to  treat 
the  node  activation  and  de-activation  is  based  on  the  notion  of  space-time  slabs. 
One  first  notes  that  the  time  integration  can  generally  be  written  in  the  form 
of  a  time  integral  that  combines  with  the  space  integration.  For  node  /,  we  will 
have  an  integral  of  the  form: 

*n+l 

/t  r  or  r 

J*j(x,t)(-j£  + ••■)<&&  (8) 

And  $/(x,£)  is  a  test  function  associated  with  node  /.  In  general,  this  test 
function  will  be  zero  where  node  I  has  no  influence  in  space  time.  Figure  5.1 
shows  for  a  one  dimensional  problem  the  active  part  of  the  space-time  slab  for 
an  element.  This  figure  also  illustrates  the  fact  that  a  node  may  be  inactive  at 
time-step  n  while  it  becomes  active  at  time-step  n+  1.  In  such  a  case,  before  we 
proceed  with  the  integration  from  tn  to  £n+1,  the  variables  associated  with  an 
inactive  node  have  no  value.  Consequently,  in  the  time  integration  between  time 
tn  and  tn+1,  for  a  previously  inactive  node,  both  U 7n  and  U?+l  are  unknown. 
For  such  nodes,  we  have  to  write  two  versions  of  equation  8  with  the  shape 
functions  and  . 

The  additional  equations,  associated  with  the  test  function  <I>y  provide  the 
required  number  of  equation  to  completely  solve  the  problem. 
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Status 


At  the  end  of  the  project,  we  had  demonstrated  that  the  node  activation/de¬ 
activation  method  worked  on  a  sample  one  dimensional  problem  and  the  imple- 
mentation  in  two  dimensions  is  still  in  progress. 
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